% % File		:	verity/ppg_0_1_display
% % Description	:   for Heart Rate Algorithm
% % Fuction		:	PPG signal display
% % Author		:	@hiyangdong
% % Version		:	V0.1
% % Time		:	12 Jan. 2016

% 0.1.1 read data and orignal calculate
clear all
orignal_txt = textread('PAH0114_1921.txt','%s');    %#ok<DTXTRD> %read txt file data 32bit with 4byte and Low byte first
orignal_dec = hex2dec(orignal_txt);                 %hex 2 dec every byte
orignal_num = length(orignal_dec);                  %length of the orignal data
% orignal_num = 1000;                  %length of the orignal data
ppg_data_byte = 20;
ppg_num  = orignal_num/ppg_data_byte;                           %length of the PPG data
ppg_data = zeros(1,ppg_num);                        %generate an 1-by-length empty array

% % 0.1.2 convert orignal hexadecimal data to ppg 32bit array
% j=1;
% for i=1:4:orignal_num
%     ppg_data(j)=(orignal_dec(i)+(orignal_dec(i+1)*256)+(orignal_dec(i+2)*256*256)+(orignal_dec(i+3)*256*256*256));
%     j=j+1;
% end

% 0.9.1 convert orignal hexadecimal data
% data format: 
% 1 byte address 
% 2 byte PPG orignal data
% 6 byte delta ppg data relative to it's fore-data --9 byte
% 1 byte Gyro X 
% 1 byte Gyro Y 
% 1 byte Gyro Z 
% 1 byte Acce X 
% 1 byte Acce Y 
% 1 byte Acce Z 
% 1 byte Magn X 
% 1 byte Magn Y 
% 1 byte Magn Z --9 byte
% 1 byte Barometer Z 
% 1 byte Temprature --2 byte

j=1;
for i=1:ppg_data_byte:orignal_num
    ppg_data(j)=orignal_dec(i+1)*256 + orignal_dec(i+2);
    j=j+1;
    ppg_data(j)=orignal_dec(i+3)*256 + orignal_dec(i+4);
    j=j+1;
    ppg_data(j)=orignal_dec(i+5)*256 + orignal_dec(i+6);
    j=j+1;
    ppg_data(j)=orignal_dec(i+7)*256 + orignal_dec(i+8);
    j=j+1;  
end

% 0.9.2 convert orignal hexadecimal data to ppg 16bit array
j=1;
m = 9;
for i=1:ppg_data_byte:orignal_num
    gyro_x_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    gyro_y_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    gyro_z_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    acce_x_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    acce_y_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    acce_z_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    magn_x_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    magn_y_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    magn_z_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    baro_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    temp_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

% % 0.4.2 ppg data max and min analysis
% ppg_data_max = max(ppg_data)
% ppg_data_min = min(ppg_data)
% ppg_data_minus = ppg_data_max - ppg_data_min
% for i=1:ppg_num-1
%     ppg_data_inc(j)=ppg_data(i+1) - ppg_data(i);
%     j=j+1;
% end
% ppg_data_inc_max = max(ppg_data_inc)
% ppg_data_inc_min = min(ppg_data_inc)
% ppg_data_inc_minus = ppg_data_inc_max - ppg_data_inc_min
% num232 = power(2,32)
% num224 = power(2,24)
% num216 = power(2,16)

% 0.1.3 time domain analysis
ppg_num = length(ppg_data);
sampletime = 0.02;          %sample time that unit is second
t = 0:sampletime:(ppg_num-1) * sampletime;      %time point

% 0.1.4 time domain plot
subplot(2,3,1);
plot(t,ppg_data,'k');
hold on
title('Time domain orignal PPG data');

% 0.1.5 frequency domain basic analysis
samplefreq=1/sampletime;    %sample frequency that unit is Hz
f=samplefreq*(0:(ppg_num-1))/ppg_num;           %frequency point

% 0.1.6 frequency domain offset calculate analysis
% sampledata_ave=0;
% for i=1:ppg_num
%     sampledata_ave = sampledata_ave+ppg_data(i);
% end
% sampledata_ave=sampledata_ave/ppg_num;
% sampledata_offset=ppg_data-sampledata_ave;

% 0.1.6.1 modify the average algorithm
sampledata_offset=ppg_data-mean(ppg_data);

% 0.1.6.2 add the plot of offset
% sampledata_offset2=sampledata_offset/30+127;
% plot(t,sampledata_offset2,'k');

% 0.1.7 frequency domain FFT analysis
fft_sampledata=fft(sampledata_offset);
fft_sampledata_abs=abs(fft_sampledata);

% 0.1.8 frequency domain plot
N=ppg_num;
subplot(2,3,3)
plot(f(1:N/2),fft_sampledata_abs(1:N/2),'k');
title('Frequency domain orignal PPG data ');

% 0.3.1 lowpass filter design
fpass = 5
fstop = samplefreq;
f_num = fix(fpass*length(t)/samplefreq);
fft_sampledata(f_num:(length(t)-f_num+1)) = 0;

% 0.3.2 lowpass filtered data
ppg_fliterdata = real(ifft(fft_sampledata)) + mean(ppg_data);

% 0.3.3 lowpass filtered data display
subplot(2,3,2);
plot(t,ppg_fliterdata);
hold on
title('Time domain filtered PPG data');

% 0.3.4 filtered data frequency domain FFT analysis
fft_fliterdata=fft(ppg_fliterdata);
fft_fliterdata_abs=abs(fft_fliterdata);

% % 0.3.5 filtered data frequency domain plot
% subplot(2,1,2);
% plot(f(1:N/2),fft_fliterdata_abs(1:N/2),'b');
% title('Frequency domain filtered PPG data');

% 0.3.5.1 filtered data frequency domain plot
subplot(2,3,4);
fft_fliterdata_abs = abs(fft_sampledata);
plot(f(1:N/2),abs(fft_sampledata(1:N/2)),'b');
title('Frequency domain filtered PPG data');

% 0.2.1 peaks and throughs detection
ds=diff(ppg_fliterdata);
filter=find(ds(2:end)==0)+1;%%find zeros
ds(filter)=ds(filter-1);%%replace zeros
ds=sign(ds);
ds=diff(ds);
throughs=find(ds>0); %when 0 change to 2 the point disappear
peaks=find(ds<0);

% 0.2.2 peaks and throughs addtional plot
subplot(2,3,2);
plot(t(peaks+1),ppg_fliterdata(peaks+1),'or');
plot(t(throughs+1),ppg_fliterdata(throughs+1),'b*');

% % 0.5.1 heart beat calculate
% ppg_time = max(t)  % t is no problem which I check with my watch
% ppg_peak = length(peaks)
% heart_rate = ppg_peak/ppg_time*60

% 0.5.1.1  heart beat calculate
ppg_peak = 0;
ppg_time = max(t)
peaks_num = length(peaks)
for k = 1:(peaks_num-1)
    if (t(peaks(k+1)) - t(peaks(k)))>0.3
        ppg_peak = ppg_peak+1;
    end  
end
ppg_peak
heart_rate = ppg_peak/ppg_time*60

% 0.8.3  heart beat calculate
subplot(2,3,5);
plot(t(1:N/4),gyro_x_data,'k');
% set(gca,'YLim',[200,300]);
subplot(2,3,6);
plot(t(1:N/4),gyro_y_data,'k');
% set(gca,'YLim',[0,10]);

